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Abstract 

We use a self-consistent Ornstein-Zernike approximation to study the Blume- 
Capel ferromagnet on three-dimensional lattices. The correlation functions 
and the thermodynamics are obtained from the solution of two coupled partial 
differential equations. The theory provides a comprehensive and accurate 
description of the phase diagram in all regions, including the wing boundaries 
in non-zero magnetic field. In particular, the coordinates of the tricritical 
point are in very good agreement with the best estimates from simulation 
or series expansion. Numerical and analytical analysis strongly suggest that 
the theory predicts a universal Ising-like critical behavior along the A-line 
and the wing critical lines, and a tricritical behavior governed by mean-field 
exponents. 

Key words: Blume-Capel model, Ornstein-Zernike approximation, tricriti- 
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I. INTRODUCTION 



The self-consistent Ornstein-Zernike approximation (SCOZA) has been introduced some 
time ago by Hoye and Stell as a method for obtaining thermodynamic and structural 
properties of simple fluid and lattice-gas systems. Like the mean-spherical approximation, 
this approach is based on the assumption that the direct correlation function C(r), which is 
related to the two-particle distribution function G{r) via the Ornstein-Zernike (OZ) equa- 
tion, is proportional to the pair potential outside the hard-core region (or, for a lattice-gas, 
for r 7^ 0) . But the dependence of the proportionality constant on density and temperature 
is determined in such a way that the same free energy is obtained from fluctuation theory 
- the so-called compressibility or susceptibility route - and from integration of the internal 
energy with respect to the inverse temperature. For the lattice-gas with nearest-neighbor 
attractive interactions (or equivalently, for the ferromagnetic spin-i Ising model), this ther- 
modynamic self-consistency is embodied in a partial differential equation whose solution, 
along with the requirement of single site occupancy, fixes C(r) (and thus G{r) ) uniquely. 
Because of numerical difficulties, this equation was only solved recently 0, showing that 
the SCOZA provides an accurate description of the properties of the 3-d Ising model over 
most of the phase diagram. The predicted values of Tc for the various cubic lattices are 
within 0.2% of their best estimates, the effective critical exponents are faithfull to the true 
behavior above except in a very narrow neighborhood of the critical point, and the zero- 
field magnetization is described asymptotically by the nonclassical exponent (3 = 0.35 [§. 
The SCOZA has been also extended to n-component and continuous spins [Q, and accurate 
results have been obtained for the hard-core Yukawa fiuid |^ and for several spin systems 
in the presence of quenched disorder [^]. 

The purpose of this paper is to apply the same type of approximation to the Blume- 
Capel model [0, a special case of the spin-1 Blume- Emery-Griffiths (BEG) model which 
represents a variety of interesting physical systems, in particular ^He -^He mixtures. This 
model has played an important role in the development of the theory of tricritical phenomena 
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TO] and has been very actively studied over the years, following the original mean field 



treatment of BEG ||^. Various methods like series expansions |]TT[, renormalizat ion-group 



calculations , and Monte-Carlo simulations |]T3[ have been used to describe the first- and 
second order regions, the tricritical region, and the crossover between them. A study of the 
coexistence curve using the mean-spherical approximation has also been proposed recently 



1^ . There is no analytical theory, however, which is able to provide a comprehensive and 



accurate description of the phase diagram in all regions (including the "wing" boundaries 
in non-zero magnetic field). As we shall see in the following, the SCOZA reaches this goal 
quite successfully whithout requiring prohibitive computational effort. In particular, the 
coordinates of the tricritical point (TCP) for the various cubic lattices are predicted with 
very good accuracy and the universal asymptotic tricritical behavior is well described. This 
suggests that the SCOZA is a reliable theory for exploring three-dimensional systems which 
exhibit first-order as well as continuous transitions. 

The paper is organized as follows: in section 2 we describe the theory and derive the 
partial differential equations that encode the thermodynamics of the model, in section 3 
we present our results for the phase diagram, and in section 4 we discuss the universal 
properties in the second-order and tricritical regions. Our conclusions are drawn in section 
5. The extension of the theory to the full spin-1 Hamiltonian is presented in Appendix A 
and details on the scaling behavior near the TCP are reported in Appendix B. 



II. THEORY 

The Blume-Capel (BC) model is defined by the Hamiltonian 

Ubc = -JY1 ^^^^ +AY,Sf-hY,S^ (1) 

<ij> i i 

where Si = 0, ±1 is the spin variable at each site i of a d-dimensional lattice and the first 
term sums over all nearest-neighbor (n.n.) pairs. This is a special case of the Blume-Emery- 
Griffiths (BEG) Hamiltonian Hbeg = T^bc ~ ^ Yli<ij> "^f ^"j i which is a microscopic 
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model for ^He He mixtures. Si = Q represents a ^He atom at site i and Si = ±1 a ^He 
atom, with the sign in the latter case describing the superfiuid degree of freedom. In this 
interpretation, the coupling constant J > is a potential that promotes superfluidity, the 
crystal-field A reflects the chemical-potential difference between the isotopes, and K repre- 
sents the difference in the van der Waals interactions between the isotopes (the actual value 
of K/J is small, so that setting K = is a sensible approximation). The magnetization 
m =< Si > identifies to the superfiuid order parameter and x = 1— < Sf > represents the 
^He concentration. As is well known, the phase diagram of ^ife — ^ He mixtures presents 
a line of second-order transitions (the so-called A-line) at high temperatures and high ^ife 
concentrations and a coexistence region associated with a first-order transition at low tem- 
peratures. The BC Hamiltonian may also describe a spin-i Ising model with a fractional 
concentration x of non-magnetic impurities in thermal equilibrium with the spin system (an- 
nealed dilution). A is then interpreted as the chemical potential that controls the impurity 
concentration (the case of quenched dilution has been studied in Ref. 0). 

Our theory for the Blume-Capel model is based on an Ornstein-Zernike approximation 
for the direct correlation function Cij which is the inverse of the connected pair correlation 
function Gtj =< SiSj > — < Si >< Sj >, i.e., 

GikCkj = 5ij (2) 

k 

where Sij is the Kronecker symbol. This OZ equation may be considered as the definition of 
Cij. It is also a consequence of the Legendre transform 

g = + ^ hiUii (3) 

i 

that defines the Gibbs free energy Q 

from the free energy JF = — A;BTlnTrexp[— 7iBc/(^B^)]- In Eq- (3), a site-dependent 
magnetic field hi has been introduced for convenience and < Si >= —dJ-'/dhi = uii is the 
local magnetization. The second functional derivatives of and Q with respect to the local 
fields and local magnetizations generate Gij and Cjj, respectively, 
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where T = PJ^,G = PG, and hi = [3hi {P = {ksT)'^ is the inverse temperature). For a 
uniform magnetic field hi = h, the system is translationally invariant and the correlation 
functions only depend on the vector r that connects the two sites. The OZ equation then 
simplifies to 

C{k)G{k) = 1 (6) 

in Fourier space. 

In contrast with G{r), the direct correlation function is expected to remain a short- 
ranged function even in the critical region (specifically, C(k = 0) < +oo). In the following, 
we shall assume that C(r) has always the same range as the pair potential. This OZ ansatz 
is the only approximation in our theory. Since the exchange interaction in the Blume-Capel 
Hamiltonian is limited to nearest-neighbor sites, this amounts to setting 

C(r) = co(J, A, m)Srfl + ci(J, A, m)5r,e (7) 

or in Fourier space, 

C{k) = co( J, A, m) [1 - z{J, A, m) A(k)] (8) 

where e is a vector from the origin to one of its nearest neighbors, A(k) = ^ Yle ^^^ '^ 
the characteristic function of the lattice (c is the coordination number), and z = — f^c. Cq 
and Ci (or, equivalently, Cq and z) are functions of J = /3J,A = f3A and m, which will 
be obtained below from the solution of the SCOZA partial differential equations (in the 
simpler mean-spherical approximation studied in ref. ||T^, one has just ci = J). It is worth 



noticing that the range of C(r) is exactly limited to n.n. separation in one dimension. 



This can be easily checked by using the transfer matrix method (see e.g. ref. [|^) to 



calculate G{r) and then inverting the OZ equation to get the direct correlation function 
(this result holds for the most general spin-1 Hamiltonian with n.n. pair interactions Ti = 
Hbc — K J2<ij> SfSj — L J2<ij> SiSj{Si + Sj) which is used as a model for ternary mixtures 
||16||). C(r) has the same spatial structure in the limit of infinite dimension (where the mean- 
field approximation becomes exact), and we expect that Eq. (8) is a reasonable assumption 
for d = 3. A major consequence is that G{r) is given in any dimension by 

G(r) = ip(r,z) (9) 
Co 



where 



g-ik.r 



P(r,z) = -J dk . (10) 

is the lattice Green's function (for instance, A(k) = ^{cosk^ + cosky + cos/c^) for the simple 
cubic lattice). The variable z (with < z < 1) is related to the second-moment correlation 
length i defined by G(k) ~ G'(0)(1 + ^^^2)^ ^ _ g, where /3G(0) = dm/dh = (3/[co{l - z)]. 
Specifically, one has 2dC,'^ = z/{l — z) for the simple hypercubic lattice. Therefore, for a 
given value of the crystal field A, the condition z = 1 gives the locus of diverging correlation 
length and diverging susceptibility in the (J, m) plane. This defines a spinodal surface in 
the (J, A, m) space; in particular, z{J, A, m = 0) = 1 corresponds to the A-line in the region 
of the h = phase diagram where the transition is continuous. 

In order to determine the two unknown functions Cq and z, we impose thermodynamic 
self-consistency. To this end, we consider the change in the free energy associated to in- 
finitesimal changes in J, A and h: 

Sjr = <: s,Sj > +6 A J2 < Sf > -6hJ2 < Si > . (11) 

<ij> i i 

In terms of the pair correlation function, this gives 

S:F/N = -^[G{r = e) + m^]6\ + [G(r = 0) + m'^]6A - m6h (12) 

where N is the number of lattice sites and the coordination number c has been adsorbed in 
the new inverse temperature variable A = cJ. The corresponding change in the Gibbs free 
energy is 



8g/N = — [Gir = e) + m^]6X + [G{r = 0) + m^]6A + Mm . (13) 
On the other hand, from Eq. (5), we have 

Therefore, in order to get the same Gibbs free energy when integrating with respect to A, A 
or m, the following Maxwell relations must be satisfied 

dC{k = 0) 1 



d\ 2 dm^ 

dG{r = 0) 1 dGir = e 



[G(r = e) + m^] (15a) 

(15b) 



dX 2 OA 

^ = —- rGr = +m2 . 15c 

dA dm^^ \ J 1 \ J 

Clearly, only two of these equations are independent and in the following we shall use Eqs. 
(15a) and (15b). 

Replacing A by the new variable r = (1 + |e^)~^ which varies from to 1, and inserting 
the explicit expressions of C(k) at k = 0, and G{r) at r = and r = e that are obtained 
from Eqs. (8-10), we finally get the two SCOZA equations 

d , , 1 Piz) - 1 , , 



dX Co 2 Ot zcq 

where P{z) = P(r = 0,2;) and the relation P(r = e,z) = {P{z) — l)/z has been used. 
Given the appropriate boundary conditions, integration of these coupled partial differential 
equations (PDE) in the two unknown functions Co(A,r, m) and z{X,T,m) gives at once the 
thermodynamics of the Blume-Capel model in the whole parameter space. Indeed, according 
to Eq. (13), the Gibbs free energy ^(A, r, m) can be obtained in the same run of integration 
from the integral of —^[G{r = e) + m^] with respect to A (thanks to the thermodynamic 
consistency, this is equivalent to integration with respect to A or to m). At fixed A and r, 
^ as a function of m has the same form as in mean- field theory except that the unstable 
region inside the spinodal is excluded. At the critical temperature, there is a single minimum 



at m — ioT T > Tt or three minima at and ±Am (with Q{±Am) — G{0)) for r < Tj. 
This defines the second- and first-order parts of the h = phase diagram, respectively. The 
maxima of the spinodal curves in the T — m plane (corresponding to d^Q /dw? = 0) define 
the lines of second-order critical points, i.e., the A-line for r > (m = 0) and the wing 
critical lines for r < rt (m = ±mc(r)). In the latter case, the corrresp ending critical field is 
given by ±hc = d{Q / N) / dm\^^j^^^. The coordinates {Tt,Tt) of the TCP can be determined 
accurately by observing the change in the convexity of the spinodal at m = 0. 

Prom a computational point of view, the coupled PDE's, Eqs. (16), define an initial 
value problem. The inverse temperature variable A plays the role of time and the equations 
describe how z{X, r, m) and co(A, r, m) propagate forward in time. We thus need to specify 
the initial condition at A = and the boundary conditions at r = 1, r = 0, and m = ±1 
(actually, because of symmetry, one can restrict the domain of integration to m > 0). 

The initial condition J — X — corresponds to the high-temperature limit where the 
spins are independent. The properties of the system can be calculated exactly and the 
correlation functions are non-zero only at r = 0. This implies that z = and C{0) = 
dh/dm — cq. The inverse susceptibility is readily obtained from the expression of the 
magnetization 

m = — — ^ ^ , (17) 

+ e'' + e-^ 

which yields 

[l-^ + (»r-2»;V + r^)V^]^ 

° (l-m2)[(l-T)(m2-2m2T + T2)V2 + (rn2_2m2r + T2)] ■ ^ 

The boundary condition at r = 1 is given by the solution of the SCOZA equation for 

the spin-| Ising model. Indeed, this limit is reached when A — oo, and the 5 = state is 

thus completely suppressed. Eq. (16b) then shows that P{z)/cq — f{m) and setting r — 1 

in Eq. (18) gives f{m) = 1 — m^. The equation for the remaining variable z{\,m) , Eq. 

(16a), becomes 

_^|_(1_,)P(„._1_1_^[(1-™=)^^1 (19) 
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which is the SCOZA equation for the Ising model studied in ref. p[ (with the standard 
replacement of lattice-gas variables by spin variables). In ref. 0, the unknown function cq 
was determined by using the hard-spin condition 5*^^ = 1 which implies that G{r = 0) = 
P{z)/cq = 1 — w? (in lattice-gas language, this is the so-called core condition). Since the 
self-consistency conditions, Eqs. (15), are exact, it it not surprising that the same result 
comes out from our equations when the state S* = is suppressed. 

The second boundary at r = is reached when A +oo. The S = ±1 states are then 
suppressed. This only happens, however, if the magnetic field h is finite. For h — > ±oo, one 
can still have a non-zero magnetization. As a consequence, z remains a non-trivial function 
of temperature and magnetization. The solution of Eqs. (16b) is again P{z)/cq = f{m), 
and setting r = in Eq. (18) gives f{m) = m(l — m) for m > 0. This leads to the equation 

Id 1 d"^ P(z) -1 

:i-^)^w = -i-i#^Hi-^)^^-^^] (20) 



m{l — m) dX 2 dm? zP{z) 

which identifies to Eq. (19) by replacing m by 1^ and A by 4A. Therefore, remarquably, 
the spinodal for h ±oo can be deduced from the spinodal of the spin-| Ising model. The 
two maxima at mc = ±| and Tc = It/'**"^ correspond to second-order transitions which 
mark the end of the wing critical lines for — > ±oo, as illustrated below. 

Finally, the boundary m = 1 is reached when all spins are in the S = 1 state. Since 
there are no more fluctuations, one has ^ = and thus 2; = 0, whereas G{y = 0) =< 5*^^ > 
— < Si >^= P{z)/cq = implies that Cq 00. 

The numerical integration of the PDE's was performed by using a finite-difference scheme 
in which the three variables A, r and m are discretized and the partial derivatives are ap- 



proximated by finite-difference representations ||T^. The first derivatives with respect to A 
are used to update z and Cq at the temperature step n + 1 by evaluating the first and second 
derivatives with respect to r and m at the step n. In principle, this is a straightforward 
procedure. There are, however, two difficulties. First, the region of integration is bounded 
by the spinodal surface which is not known in advance. Secondly, there is a singular behav- 
ior as one approaches the spinodal. To see this, let us rewrite the PDE's using as unknown 
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functions z and the new variable v = P{z)/cq. We get 

|^ = ^^(l-r)|:(^Wt;) (21b) 

where A{z,v) = -Piz){l - z)/v^, B{z,v) = il/v)d[{l - z)P{z)]/dz and i){z) = {P{z) - 
1) / {zP{z)). Eq. (21a) can be viewed as a nonhnear diffusion equation and plays 
the role of a diffusion coefficient that diverges like (1 — z)~^ when z ^ 1, namely on the 
spinodal surface. These two difficulties are already present in the equation for the Ising 
model, Eq.(19), but then the spinodal is just a line in the (A, m) plane 0. 

We solved the first problem as follows. Whenever the variable 2;„(t, m) at the temperature 
step n enters the interval (1 — e, 1), where e is small number (typically, e < 10^^ — 10~^), 
we consider that the spinodal is reached and we stop the calculation. The spinodal is then 
defined by the corresponding values of r and m. At the next temperature step, we use 
the same values Zn and f„ to compute the partial derivatives with respect to r and m on 
the spinodal (in other words, we locally "freeze" the values of z and v inside the spinodal). 
On the other hand, the vanishingly small values of A{z,v) for z 1 impose a dramatic 
decrease of the inverse temperature spacing AA as the spinodal is approached. Indeed, as 



is well known solving a diffusion equation by an explicit method requires to keep the 
"time" step below a certain value which is proportional to D^^, where D is the diffusion 
coefficient. In ref. 0, this problem was avoided by using an implicit method which is 
inconditionally stable. Unfortunately, it is not easy to generalize this procedure to a system 
of nonlinear PDE's like Eqs. (21) and we had to keep a simple explicit algorithm. In the high- 
temperature region, we typically adopted the spacings Am = At = 10^^ and AA = 10^'^. 
In the vicinity of the A-line and in the tricritical region, Ar was set at 2.10^^ and AA was 
gradually decreased down to 10"''. This allowed to determine the critical parameters with 
excellent accuracy. For instance, in the limit r — > 1, we found Jc = 0.22125 for the inverse 
critical temperature of the Ising model on the simple cubic lattice. This corresponds to 
J^'^ = AJc = 0.88500 for the n.n. lattice-gas, in perfect agreement with the value obtained 
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in ref. p[ (this is also within 0.2% of the best-estimate result [0). When higher accuracy 
was required, for instance to determine the asymptotic critical exponents, AA was further 
decreased to 10~^. The integration was usually carried down to fc^T/Jc ~ 0.18 (below this 
value, the spinodal lines in the vicinity of r = 1 are so close to the boundary m = 1 that 
we could not integrate the equations with sufficient accuracy while keeping a reasonable 
spacing Am; indeed, decreasing Am forces to decrease also AA in order to avoid numerical 
instabilities [|17| ). 

Before presenting the numerical results, it is instructive to consider the high-temperature 
series expansion of the solution and compare it to the exact results. Since z — > for 
A ^ 0, one can replace the Green's function P{z) by its expansion in powers of z. We then 
expand z and Cq as double series in A and m, z{\, r, m) = J2pq ^pqij)'^^'^'^'^ ^"^^ Co(A, r, m) = 
Ylipq c^q{j)\^rn?'i . The coefficients Zpqij) and c^qij) are polynomials of r that satisfy a system 
of linear algebraic equations at each order in A and m. In the case of the fee lattice (c = 12, 
P{z) = 1 + z^/12 + + 52;'^/192 + 52:^/288 + ...) for which extensive series expansions 

have been derived by Saul et al. |]TT|, the results for the zero-field ordering susceptibility 



Xo = Xlr ^(^' = 0) = dm/ dh\j^^Q and the second moment of the correlation function 
= Er ^^^(r; m = 0) = c^^Xo are 

kBTxo = r +12t^J + 6{t^ + 21t^)J^ + 2{r'^ + l^r" + 621r^)P + ^(^^ + 234r3 + 5115r^ 
+ 23778r^) + ^^(r^ + 612r^ + 31851r^ + 342690r^ + 1122462r^) 
+ . . . (22) 

and 

/i2 = 12r2 J + 288r3p + 2{t^ + 66r^ + 2385r^)P + ^(240r^ + 10080r^ + 133488r^) 

+ ^(r^ + 492r3 + 50931r^ + 1156410r^ + 8474742r^) + . . . (23) 

Comparison with the exact series expansions shows that the coefficients in both expressions 
are exact through order {l3jy while higher-order terms do not deviate significantly from the 
exact ones. This is similar to the case of the spin-i Ising model and we take this result as a 
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strong indication that the numerical predictions of the SCOZA should be very close to the 
exact solution of the model. 

To close this section, let us note that the present theory does not provide a complete 
description of the system. In particular, it does not give any information concerning the 
two correlation functions Gfj^^ =< S^Sj > - < Sf >< Sj > and Gff^ =< SiSj > - < 
Si >< Sj >. In order to determine these functions, one needs to introduce a set of three 
different direct correlation functions. This implies to perform a double Legendre transform 
that defines a new Gibbs free energy which is a function of m and x, the concentration order 
parameter, instead of m and A. The main interest of this alternative theory is that it allows 
to study the general spin-1 Hamiltonian with K ^ and L ^ 0. On the other hand, the 
numerical solution is more difficult as one has to solve three coupled PDE's instead of two. 
Details on the derivation of these equations are given in Appendix A. 



III. RESULTS 

In this section we concentrate on the SCOZA numerical predictions for the phase bound- 
aries. These are nonuniversal properties which are lattice-dependent. If not stated otherwise, 
the results presented here correspond to the simple cubic lattice for which no systematic 
study has been performed in the literature. 

The overall shape of the spinodal surface in the (T, r, m) space and the vicinity of the 
TCP are depicted in Figs. 1 and 2, respectively. We clearly see the evolution from the single 
curve at r = 1 (the spinodal of the spin-| model) which has a maximum at m = to the 
two symmetrical curves at r = with maxima located at = ±|, marking the end of the 
wing critical lines. 

The Tc(r), Tc(A), and Tc{x) phase diagrams in zero field are shown in Figs. 3-5. Second- 
and first-order phase boundaries are shown as full and dashed lines, respectively. The curves 
are quite similar to those obtained by series expansions [|TT]| and Monte-Carlo simulations 



for the fee lattice. In particular, we see that the slope of the phase boundary across the TCP 
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is finite and continuous in both Tc(r) and T'c(A) (specifically, we find (Tt/Jc) dA/dT\j.^ = 
—0.045); the Tc(A) phase boundary is slightly concave upward just below the TCP, and the 
A-line appears to extrapolate into the interior of the two-phase region in the Tc{x) phase 
diagram (a continuous slope, however, cannot be strictly ruled out by our calculations). As 
is well known, the slope of the A-line and the slope of the coexistence curve on the ^He- 
rich side are not the same experimentally. This is also predicted by renormalizat ion-group 
analysis f^, in contrast with mean-field theory P]. 

The accuracy of our calculation for the A-line in the A — T plane can be checked for 
the special value A = ln2 for which a careful Monte-Carlo calculation and finite-size study 



has been performed by Blote et al. ||T8|. Our prediction for the inverse critical temperature 
Jc = J/lksTc) = 0.3924 is in excellent agreement their estimate = 0.3934224(10). The 
accuracy of the theory is thus the same as for the spin-i Ising model p|. 

As noted earlier, the TCP corresponds to the point where the convexity of the spinodal 
in the T — r plane changes at m = 0. This is illustrated in Fig. 6 which shows the 
temperature dependence of the order parameter and the spinodal lines in the tricritical 
region on the first-order side of the phase boundary (observe that Am(r), the discontinuity 
in the order parameter across the first-order phase boundary, moves away from the spinodal 
as r decreases). The coordinates of the tricritical point are ksTt/ J = 1.4160 ± 0.0040, 
Tt = 0.2114 ± 0.0010 {At/ J = 2.8457), Xt = 0.655 ± 0.006, where the uncertainties reflect 
the flnite size of the grid spacings. The predictions for Tf and A^ are in excellent agreement 



with the recent Monte-Carlo estimates of Deserno |T3|: ksTf/J = 1.4182 ± 0.0055, At/ J = 
2.8448 ± 0.0003 (these numbers, however, are different from those quoted in Ref. [|19| which 
locate the TCP near ksTt/J = 1.3900, At/ J = 2.849, and Xt = 0.61; if these (unpubhshed) 
results are correct, our value of Xt is overestimated and too close to the mean-fleld prediction, 
xf ^ = 2/3). Similarly, for the fee lattice, we flnd keTt/J = 3.1116 ± 0.0090, n = 0.2454 ± 
0.0010 {At/ J = 5.6520), xt = 0.658 ± 0.006, which we may compare with the Monte- 
Carlo estimates of Jain and Landau keTt/J = 3.072 ± 0.024, At/ J = 5.652 ± 0.048, 
Xt = 0.56 ± 0.02 (note that our value of Xt is in much better agreement with the series 
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expansion estimate of Saul et al. |TT|, xt = 0.6651o;q°5; obviously, further work is needed to 
locate precisely the TCP in the x — T plane). Finally, our predictions for the bcc lattice are 
kBTt/J = 2.0264 ±0.0060, n = 0.2354 ± 0.0010 {At/ J = 3.7918), Xt = 0.656 ±0.006 (to our 
knowledge, this lattice has only been studied by real-space renormalization-group methods 



which do not predict accurately the location of the tricritical point). 
At the TCP, the A-line bifurcates into two symetrical wing critical lines. The projections 
of the wing boundaries onto the A — T, A — h, and T — h planes are shown in Fig. 7. 
Mean-field theory |^ predicts that the critical field should go to infinity at ksTc/Jc = |. 
We clearly see in Fig. 7 that this value is overestimated. In fact, as noted earlier, the present 
theory predicts that he ±oo at kBTc/ Jc = |/cbT/**"^/ Jc = 0.188. For the fee lattice, this 
yields ksTc/Jc = 0.204, which is consistent with the value that can be extracted from the 



Monte-Carlo simulations of Jain and Landau 13 



IV. ASYMPTOTIC BEHAVIOR IN THE CRITICAL AND TRICRITICAL 

REGIONS 

As mentioned in the Introduction, the SCOZA for the 3-d spin-| Ising model has a 
nontrivial scaling behavior in the critical region [0,0]. Above Tc, the asymptotic behavior is 
the same as in the mean-spherical approximation and the exponents are those of the spherical 
model. This spherical scaling, however, is detectable only in a very narrow neighborhood 
of the critical point, and the effective SCOZA exponents are close to the true Ising ones 
down to reduced temperatures of around 10~^. Below Tc, the scaling is neither spherical nor 
classical with two scaling functions instead of one 0. Despite this shortcoming, the zero- 
field magnetization is very well described, with an asymptotic exponent /5 = 7/20 = 0.35 
which is close to the exact value (3 ~ 0.33. It is thus interesting to also investigate the 
critical behavior of our SCOZA equations in the critical and tricritical regions of the 3-d BC 
model. 

We first consider the behavior of the zero- field ordering susceptibility xo as T ^ ^c(t) 
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along paths of constant r for r > Tj. Accurate evaluations were relatively straightforward 
to perform in the disordered phase: we only had to gradually decrease the spacing A A 
as discussed earlier. Fig. 8 shows a log-log plot of ksTxa as a function of the reduced 
temperature t — 1 — y' together with the corresponding effective exponent 7e// defined as 
the local slope d \og{k bT xo) /d log t. In the region 1.0 > r > 0.25, it can be seen that each 
7e//(T) reaches the value 2 for t ~ 10~^ as in the case of the spin-| model, showing that the 
asymptotic spherical behavior is universal. This is no more true when one moves further 
away from Tc. However, in the range t > 10~^, a quasi-universal behavior is still observed 
for 1.0 > r > 0.6, and the critical behavior is governed by an effective exponent which is 
close to the exact Ising value 7 1.24. On the other hand, as r approaches its tricritical 
value Tt — 0.211, there is an abrupt crossover to another behavior which is governed by the 
exponent jeff ~ 1 over a wide range of temperatures. There is good numerical evidence 
that Jeff reaches 1 asymptotically at r = r^. 

For subcritical temperatures, it was more difficult to obtain accurate results in the vicin- 
ity of Tc because of our use of an explicit method to integrate the PDE's. Accordingly, 
we were only able to explore the critical behavior in a restricted range of temperatures 
t — 1 — ^. Log- log plots of the temperature dependence of the order parameter m are 
shown in Figs. 9 and 10. Despite the limited range, it appears from Fig. 9 that in the 
second-order region well above the crossover to tricritical behavior, the slope of each curve 
has a common asymptotic limit which corresponds to the SCOZA prediction for the Ising 
model, Picma = 7/20. Fig. 10 shows that for smaller values of r, one needs to go closer to 
to reach this asymptotic universal regime. Again, we observe an abrupt crossover to another 
behavior as one enters the tricritical region. Our results are consistent with the asymptotic 
exponent /^^ = 1/4 for r = r^. 

Finally, we analyze the shape of the wing critical boundary as it approaches the TCP. 
Figs. 11 shows the log-log plots of the critical field he and the magnetisation mc as a function 
of the reduced temperature 1 — ^- Our numerical data arc consistent with asymptotic power- 
law behaviors governed by the exponents 5/2 for he and 1/2 for nie- 
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All the above numerical results strongly suggest that our theory describes the whole 
critical portion of the phase boundary by the same exponents as the SCOZA equation for 
the spin- 1 Ising model and that near the TCP there is a crossover to a tricritical behavior 
described by mean-field exponents. This is supported, and can be further rationalized, by 
considering a heuristic scaling analysis of the coupled SCOZA PDE's, Eqs. (16). The 
argument is summarized below and detailed in Appendix B. 

Let us assume that close to the TCP, the singular part of the Gibbs free energy can be 
written as 

Qs^n, ^ ^) (24a) 



(24b) 

where we have introduced the two scaling fields t = (T — Tt)/Tt and g = {t — Tt)/Tt — at 
where {Tt/Tt)a = dTx/dT\j,^ > is the slope of the A-line at the TCP (this is also the slope 
of the triple line below Tt, as we have seen that the slope of the phase boundary is finite 
and continuous at the TCP). Eqs. (24) have the form of the standard tricritical scaling 
hypothesis except that we use the magnetization m instead of the magnetic field h 



as variable. and Qj^^-^ are the scaling functions, where the subscript (±) represents the 
sign of t, i.e., denotes when the temperature is above or below the tricritical temperature T^. 
When |t| ^ with g = 0, the TCP is approached tangentially to the phase boundary in the 
symmetry plane h = (and Eq. (24a) is then the convenient form of the scaling hypothesis), 
whereas the TCP is approached with a finite angle with the phase boundary when g ^ 
(and Eq. (24b) is the convenient scaling form). This defines two sets of exponents (a,0, /?) 
and {at, 4>t, Pt) that are related through 2 — at = (2 — a)/0, (pt = l/4> and Pt = l3/(p- 

Since SCOZA is thermo dynamically self-consistent, the scaling behavior of the Gibbs free 
energy near the TCP is inherited by its various derivatives. In particular, the asymptotic 
behavior of the magnetic field h = d{Q / N) / dm is 

h ^ W-''-^^Q±{u,v) ^ \g\^-^-P^l-g^^{ut,Vt) (25) 
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where u = g/\t\'^, v = m/\t\^ {ut = t/\g\'^'^, Vt = m/\g\^'^), and the inverse ordering suscepti- 
bihty Xo^ = d"^ {Q / N) / dm? and the singular part of the ^He concentration order parameter 
p = 1 — X = d{Q /N)/dA obey (up to irrelevant multiplying factors) 

Xo' - ItV^S^M ^ \9\''^^gi{u,,v,) (26) 

and 

Psing ^ -\tr--'^^g^iu,v) ^ -Igr'-^-^^^Giiu^v,) (27) 

where 7 = 2 — a — 2/3 (resp. 7f = 2 — at — 2/3*). Note also that because m is used in 
Eqs. (24) instead of h, the zero-field magnetization is solution of the implicit equation 
h = d{Q/N)/dm = 0. This implies that the singular part of this quantity near the TCP 
obeys m,,„, ^ \tfM±ig/\t\*) ^ \gf^Mi{g/\t\*^). 

From the numerical results shown in Figs. (8-11), it appears that 7t ~ 1 (see the 
curve r = 0.22 in the lower part of Fig. 8), jSt ~ 1/4 (see the curve r = 0.21 in Fig. 10), 
2 — a — (3 ~ 5/2 and /5 ~ 1/2 (upper and lower parts of Fig. 11). (We have also good 
numerical evidence that /3 ~ 1/2 from a log-log plot of the discontinuity of the zero-field 
magnetization as a function of (Tj — T)/Tj across the first-order phase boundary.) All these 
exponents have their classical values, and from the scaling relation a + 2f3 + 'y = 2 (resp. 
at + 2(3t -|- 7t = 2), we deduce that 7 = 2 and at = \. These values are all consistent with a 
crossover exponent = 2. 

We show in Appendix B that the scaling ansatz, Eqs. (24) or Eqs. (26-27), is compatible 
with the asymptotic behavior of the PDE's, Eqs (16), in the tricritical region. This analysis 
indicates that a non-trivial scaling is found when the exponents obey the two relations 
7 = and 7 = 4/3, which are satisfied by the classical values. Morever, one finds that the 
scaling function of the zero-field susceptibility above Tt obeys an equation which is quite 
similar to the asymptotic SCOZA equation studied in Ref. [Q for the spin-| model. It can 
be inferred that the critical behavior of the present theory along the A-line is the same as 
the SCOZA prediction for the Ising model. This is consistent with the exponent f3 = 7/20 
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which is observed numerically in Fig. 9 along the high-temperature part of the A— line. This 
universality appears to hold along the wing critical lines too since the boundary condition 
to Eqs. (16) at the end of these lines (for he ±00) is again the SCOZA PDE for the Ising 
model, as explained in section 2. 



V. CONCLUSION 

The present study shows that a thermodynamically self-consistent OZ approximation 
provides a very good description of the properties of the 3-d Blume-Capel model in all parts 
of the phase diagram. Like in the case of the Ising model, non-universal properties such as 
the shape of the phase boundaries and the location of the tricritical point are predicted with 
remarkable accuracy. Moreover, there is good numerical and analytical evidence that the 
SCOZA correctly predicts a universal critical behavior along the A-line and the wing critical 
lines (with a zero-field magnetization exponent 0.35 that is very close to the true Ising value), 
as well as a crossover to tricritical behavior governed by classical exponents. Therefore, the 
SCOZA proves to be a powerful tool for studying spin systems which exhibit first-order 
and/or continuous transitions. This is confirmed by further work on the ferromagnetic 
spin-3/2 and Potts [H models. 
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APPENDIX A: 



In this Appendix, we derive the SCOZA equations for the the most general spin-1 Hamil- 
tonian with n.n. couphngs, 

n = nBc-KY,sts^^-LY, SiS,{s, + s,) , (ai) 

<ij> <ij> 



which is a model for ternary mixtures ||T6[. The solution of these equations also provides 
a complete description of the pair correlation functions of the Blume-Capel model. These 
functions can be generated by introducing site-dependent fields hi and Aj in the Hamiltonian 
(Al), which yields 



Gff = < S,S, >-<S,>< S, >= (A2a) 



Gff = < SiS^ >~<Si>< >= (A2b) 



Gff' = < SfS^ >-<S^>< S'' >= £^ . (A2c) 



We then perform a double Legendre transform that takes the fields hi and Aj into rrii 
and Xi, respectively, where Xi = 1— < Sf >= 1 — dJ-'/dAi. This defines a Gibbs free 
energy Q{T, {nii}, {xi}) = JF + hiirii — A,j(l — Xi) which satisfies hi = dQ/drrii and 
Aj = dQ/dxi. Q is the generating functional of the direct correlation functions. 



ss 

drriidm ^ "^^^ 



- (A3b) 

that are related to the G"s via a set of Ornstein-Zernike equations. In the limit of uniform 
fields, these equations in Fourier space take the form C(k)G(k) = 1, where G(k) and C(k) 
are symmetrical square matrices. This readily yields 



G^fk) = '-^ 2 (A4a) 
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G^^(k) = 2 (A4b) 

G^^\k) = (A4c) 

We now assume that the range of the direct correlation functions is hmited to n.n. 
separation, i.e., 

^^^(k) = - CssXm (A5a) 
(7^^'(k) = Co^^^[l-C552A(k)] (A5b) 
C'''\k) = cf'\l - Cs^s^X{k)] (A5c) 

where the cq's and the ^'s are functions of T, m and x to be determined. This fixes the form 
of the correlation functions in r-space, and after some calculations we find 

where P{z,y) is the lattice Green's function defined by Eq. (10), and ^1,^2 are related to 
the co's and the ^'s via the relations 

Css + Cs^s^ = 2R\ss^ + {1-R^) (zi + Z2) (A7a) 
Css C5252 = RVss^ + (1 - R'')ziZ2 (A7b) 

and 

it! _ [z,P{z,) - z^Pjz^) - CssjPjzi) - P{z2))mz,P{z,) - z^Pjz^) - CsMPjzi) - P{Z2))V/' 

Ro ZiP{zi) - Z2P{Z2) - CSS<P{Z1) - P{Z2)) 

(A8) 

where R — Cq^^ /{cq^Cq^^^Y^'^ and Rq is the high-temperature limit of R which can be 
calculated exactly, as explained below. Note that in Eqs. (A6) we have eliminated the cq's 
to introduce the on-site values of the correlation functions which are simple functions of 
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the order parameters m and x. Indeed, since Si can only take the values 0, ±1, one has 

< Sf >=< Si > and < Sf >=< Sf >,so that 

^^^(r = 0) = 1 - X - (A9a) 

^■^^'(r = 0) = ma; (A9b) 
G^'^\r^O)^x{l-x) . (A9c) 

In terms of these variables, one also has 

^ ' (l-i?2)(l-a;-m2) ^1-^2 

^ ^ (1 - R'^)xm zi- Z2 ^ ^ 



It remains three unknown functions to be determined, and it is convenient to choose 
Zi,Z2 and R and to use Eqs. (A7) and (A8) to calculate the C's- The three additional 
equations that we need are obtained by imposing thermodynamic self-consistency. On the 
one hand, the enthalpy is given by 

^ = -^[G'Hr = e) + m'] - ^[G'''\v = e) + (1 - x)'] - Lc[G''\v = e) + m(l - x 



On the other hand, we have from Eqs. (A3) 



(All) 



C'^k = 0) = 1^ (A12a) 



C^'"(k=0) = -^ (A12b) 



dmdx 
dx"^ 



C^'^'ik = 0) = g . (A12C) 



This yields the three Maxwell equations 





0) 


dX 






0) 


dX 






0) 


dX 



1 ^^[^^^(r = e) + QiC;^'^'(r = e) + 2a2G'^^\r = e)] 

ld^[G^^(T = e)+aiG^'^\r = e) + 2a2G^^\r = e)] 

= -"^ + 2 d^x 

1 ^^[^^^(r = e) + Q;iG^'^'(r = e) + 2a2G^^"{r = e)] , , ^ ^ , 
= -«i - 2 (Al3c) 
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where A = (5Jc, a-i — Kj J, and 0:2 = Lj J. 

These equations, together with Eqs. (A7) and (A8), constitute a set of three PDE's in the 
unknown functions Zi{\m,x), Z2{X,m,x), and R{X,m,x), whose solution encodes the full 
thermodynamics of the model Hamiltonian (Al). The initial conditions for J = X = L = 
are easily obtained since the correlation functions are then non-zero at r = only. One has 
zi = Z2 = and from Eqs. (A4) and (A9) 

mx 

° " ~[x(l-x)(l-x-m2)]V2 ■ ^^^^^ 

It sould be noticed that in the case of the Blume-Capel model (for which ai = q;2 = 0), 
this SCOZA is different from the one presented in the main text. This can checked for 
instance by computing the high-temperature expansion of the solution. Both theories yield 
zero-field properties which are exact through order A'*. It is unclear which one provides the 
best numerical predictions. 
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APPENDIX B: 



In this appendix, we show that the tricritical scahng hypothesis, Eqs. (24) or (26-27), 
is consistent with the asymptotic behavior of the SCOZA PDE's, Eqs. (16). The notations 
are those of the main text. 

To this end, it is convenient to rewrite the PDE's in terms of the two variables £ = 
(1 — z)^/^ and p = 1 — X. £"^18 proportional to the inverse susceptiblity Xo^ and p = G{v = 

0) + = P{z)/cq + m^. At the TCP, we have £ = and p = pt- For t — 0, m — > 0, and 
p Pt, Eqs. (16) take the asymptotic form 

PtXt dt ^ 2 ^^P{l)-ldm^ dm^^ ^ ' 

where t = {T — Tt)/Tt, Xt = cJ/^ksTt), 6t = {T — Tt)/Tt, and we have used the expansion of 
the 3-d lattice Green function for z — > 1, P{z) ~ -P(l)[l — b£ + 0{£'^)], where 6 is a positive 
constant In these equations and in the following, all derivatives are taken at the TCP. 
By suitably rescaling the variables as £^ ^ bpt/{P{l) — 1) £, t ^ Xtb'^pf/[P{1){P{1) — 

1) ]2 t,m^ m/(P(l) - 1)1/2 - 1)^(1 - n)] St, we obtain the two 
simplified equations 

dt 2 dw? ^ ' 

dp _ 1 d{£ - p) 

m - 2^6^ ^^^^^ 

We now introduce the tricritical scaling ansatz for £ and for the singular part of p, 
according to Eqs. (26) and (27), 

£ ^\t\^/^E±{u,v) (B3) 

Psing -\t\^-''-'^R±iu,v) (B4) 

where u = and v = Because of thermodynamic self-consistency, E± and R± obey 
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9El d'R± 

For the sake of simplicity, we only consider the case t > (i.e., T > Tt), but a similar 
analysis can be performed for t < 0. Then Eqs. (B2) yield 

(B6a) 

where we have used the scaling relation — 2 — a — 2/3 which results from thermodynamic 
self-consistency. If the cross-over exponent is greater than 1 (which is usually the case 
and is indeed found numerically), the first term in the left-hand side of Eq. (B6a) may be 
neglected asymptotically. A non-trivial scaling is then found when the exponents are related 
through the two relations 7 = = 4(3. For the same reason we may neglect the first term 
in the left-hand side of Eq. (B6b) and we obtain the relation 7 = 2(2 — a — 0) (which is not 
independent from the preceding ones). Actually, we expect that as in the SCOZA for the 
Ising model, the enthalpy is analytic in and T — when approaching a critical point 
from a disordered phase, which corresponds to 7 = 2 and (3 — 1/2. (At the tricritical point, 
the term of course vanishes.) The scahng functions satisfy the two non-trivial PDE's 



dEl _^ l d\E+- R+) 

du 2 dv 

dE+ _ , ^^^dR. 



du [l + 2a(0-l)]^ (B7b) 
Using Eq. (B5) to eliminate one of the functions, we finally obtain a single equation for £'+ 

^^^"^^ ^^^^ 
where we have used the rescahng u — > 2/[2a(0 — 1) -|- 1] u. (Note that the multiplying 
factor is positive since a > and > 1.) Of course, this equation must be accompanied by 
some boundary conditions. These are obtained from the analytical requirements that the 
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scaling functions must satisfy near the TCP and in the vicinity of the critical lines and the 



coexistence surfaces (see, e.g., Ref. |T^). It can be shown that these boundary conditions 
are also compatible with the SCOZA equations. The other function can be obtained 
as [1 + 2a{(f) — t>) = E+{u,v) — + R+o, where i?+o is a constant that can be 

determined by the boundary conditions. 

Eq. (B8) has an important consequence for the scaling behavior near the TCP when one 
approaches the A-line along a path at fixed r. The convenient temperature variable is then 
t = [T — Tx{T)]/Tt which measures the distance from A-line at fixed r, so that t = t + tx, 
where = (^a(''") — Tt)/Tt defines the A-line near the TCP. When \i\ — > 0, or t — > t\{T), the 
scaling field g behaves as g — go i, where go = {t — Tt)/Tt is a constant, and the scaling 
variables u and v behave as -u — mq ~ t and w ~ m, where uq is a constant. As a consequence, 
Eq. (B8) can be rewritten as 

This is precisely the asymptotic SCOZA equation for the spin-| Ising model which has been 
studied in Ref. (with E'^ playing the role of the variable (p in that reference). Therefore, 
one expects that the critical behavior above and below the critical temperature Tx{t) will 
be identical to that of the SCOZA for the Ising model (with for instance /? = 7/20 for the 
zero-field magnetization m{t)). 
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FIGURE CAPTIONS: 



Fig. 1: SCOZA spinodal surface of the 3-d Blume-Capel in the {T — t — m) space. 
Fig. 2 : Detail of the spinodal surface near the TCP. 

Fig. 3: Zero-field T'c(t) phase diagram. Second-order and first-order parts of the phase 
boundary are shown as full and dashed lines, respectively. For numerical reasons, calcula- 
tions have not been performed below A;sT/(Jc) ~ 0.18 (see text). 

Fig. 4: Zero-field Tc(A) phase diagram. Second-order and first-order parts of the phase 
boundary are shown as full and dashed lines, respectively. The inset describes the vicinity 
of the TCP. 

Fig. 5: Zero-field Tc{x) phase diagram. Second-order and first-order parts of the phase 
boundary are shown as full and dashed lines, respectively. The inset describes the vicinity 
of the TCP. 

Fig. 6: Temperature dependence of the order parameter (full lines) and spinodal lines 
(dotted) in the vicinity of the TCP on the first-order part of the phase boundary. 

Fig. 7: Projections of the wing boundaries onto the A — T, A — /i and T — h planes. 

Fig. 8: Log-log plot of the zero-field ordering susceptibility ksTxo as a function of the 
reduced temperature 1 — Tc/T and corresponding effective exponent 7e//- The different 
curves correspond to r = 1.00, 0.79, 0.60, 0.40, 0.29, 0.25 and 0.22. 

Fig. 9: Log-log plot of the order parameter m as a function of the reduced temperature 
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1 — T/Tc in the second-order region well above the crossover to tricritical behavior. 

Fig. 10: Log-log plot of the order parameter m as a function of the reduced temperature 
1 — T/Tc in the second-order and tricritical regions. 

Fig. 11 : Behavior of the wing boundaries as the TCP is approached: log- log plots of the 
critical field he and magnetization rric as a function of the reduced temperature 1 — Tc/Tt. 
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